MODULE pfmm_Vff_driver

!   USE fmm_stats
   USE fmm_global_consts
   USE fmm_global_types

   USE pfmm_periodic_images, ONLY: pfmm_get_periodic_image

   USE fmm_T_contractors,  ONLY: fmm_set_T_con_ptrs,               &
                                 fmm_select_T_con,                 &
                                 fmm_init_T_contractors,           &
                                 fmm_free_T_contractors
   USE fmm_W_contractors,  ONLY: fmm_init_W_contractors,           &
                                 fmm_free_W_contractors
   USE fmm_T_pair_builder, ONLY: fmm_init_T_pair_builder,          &
                                 fmm_close_T_pair_builder,         &
                                 fmm_gen_local_T_pairs,            &
                                 fmm_gen_nonlocal_T_pairs

!   USE fmm_box_utils,      ONLY: fmm_deepest_level
   USE pfmm_box_builder,   ONLY: pfmm_get_box_qlm_at_level,         &
                                 pfmm_get_box_paras_at_level
   USE fmm_W_pair_builder, ONLY: fmm_get_raw_Vff_from_boxed_Vff

   IMPLICIT NONE
   PRIVATE
   ! Public procedures
   PUBLIC :: pfmm_get_Vff

   ! Levels in box hierarchy 
   ! 1=> 9**3 supercells; 2=> 27**3 supercells
   INTEGER(INTK), PARAMETER :: pfmm_max_level = 12

CONTAINS

!-------------------------------------------------------------------------------

   SUBROUTINE pfmm_verify_data_in(RHS_in)

      IMPLICIT NONE
      TYPE(raw_mm_data), INTENT(IN) :: RHS_in
      LOGICAL :: fail

      fail = .FALSE.
      IF (.NOT.ASSOCIATED(RHS_in%paras)) fail = .TRUE.
      IF (.NOT.ASSOCIATED(RHS_in%qlm_T)) fail = .TRUE.
      IF (.NOT.ASSOCIATED(RHS_in%qlm_W)) fail = .TRUE.
      IF (SIZE(RHS_in%paras) /= SIZE(RHS_in%qlm_T,2)) fail = .TRUE.
      IF (SIZE(RHS_in%paras) /= SIZE(RHS_in%qlm_W,2)) fail = .TRUE.
      IF (fail) CALL fmm_quit('mms pointers sent incorrectly to fmm_Vff_driver')

   END SUBROUTINE pfmm_verify_data_in

!-------------------------------------------------------------------------------

   SUBROUTINE pfmm_get_Vff(scheme,cell_size,LHS_paras,RHS_mms,Vff)

      USE pfmm_box_builder, ONLY: pfmm_init_box_builder, pfmm_free_box_builder 

      IMPLICIT NONE
      TYPE(fmm_scheme),    INTENT(INOUT) :: scheme
      REAL(REALK),         INTENT(IN)    :: cell_size(3)
      TYPE(raw_mm_paras),  INTENT(INOUT) :: LHS_paras(:)
      TYPE(raw_mm_data),   INTENT(INOUT) :: RHS_mms
      REAL(REALK), TARGET, INTENT(OUT)   :: Vff(:,:)

      WRITE(IO1,*) "Computing multipole potential..."
      CALL pfmm_verify_data_in(RHS_mms)
!FIXME pass coordinate centre of supercell as an argument here
!      IF ( DOT_PRODUCT(cntr, %cntr) > 0 ) &
!         CALL fmm_quit('pfmm centre of original supercell!')
      CALL pfmm_init_box_builder(LHS_paras,RHS_mms,pfmm_max_level,cell_size)

print *, "trans_LMAX", scheme%trans_LMAX
print *, "raw_LMAX", scheme%raw_LMAX
!print *, "SIZE(RHS_mms%qlm_T,1)", SIZE(RHS_mms%qlm_T,1)
!print *, "SIZE(RHS_mms%qlm_T,2)", SIZE(RHS_mms%qlm_T,2)
      scheme%phase = FAR_FIELD
!      SELECT CASE (scheme%algorithm)
!         CASE (DO_NULL)
!            CONTINUE
!         CASE (DO_FQ)
!            print *, "Doing FQ algorithm"
!            CALL pfmm_get_FQ_Vff(scheme,cell_size,LHS_paras,RHS_mms,Vff)
!         CASE (DO_BQ)
!            print *, "Doing BQ algorithm"
!            CALL pfmm_get_BQ_Vff(scheme,cell_size,LHS_paras,Vff)
!         CASE (DO_FMM)
            print *, "Doing FMM algorithm"
            CALL pfmm_get_FMM_Vff(scheme,cell_size,LHS_paras,Vff)
!         CASE DEFAULT
!            CALL fmm_quit('invalid algorithm requested!')
!      END SELECT

      CALL pfmm_free_box_builder

   END SUBROUTINE pfmm_get_Vff

!-------------------------------------------------------------------------------

   SUBROUTINE pfmm_get_FQ_Vff(scheme,cell_size,LHS_paras,RHS_mms,Vff)

      IMPLICIT NONE
      TYPE(fmm_scheme),           INTENT(IN)    :: scheme
      REAL(REALK),                INTENT(IN)    :: cell_size(3)
      TYPE(raw_mm_paras), TARGET, INTENT(IN)    :: LHS_paras(:)
      TYPE(raw_mm_data),          INTENT(IN)    :: RHS_mms
      REAL(REALK),        TARGET, INTENT(INOUT) :: Vff(:,:)

      INTEGER(INTK), PARAMETER :: pair_type = LHS_raw_RHS_raw
      TYPE(raw_mm_paras), ALLOCATABLE :: ref_paras(:)
      TYPE(gen_mm_paras) :: LHS, RHS
      TYPE(fmm_range) :: ncell
      REAL(REALK) :: T0, T1, cntr(3)
      INTEGER(INTK) :: i,j,k,p

      CALL CPU_TIME(T0)
      NULLIFY(LHS%raw_paras,LHS%box_paras)
      NULLIFY(RHS%raw_paras,RHS%box_paras)
      LHS%raw_paras => LHS_paras(:)
      RHS%raw_paras => RHS_mms%paras(:)
      ALLOCATE(ref_paras(SIZE(RHS_mms%paras)))
      ref_paras(:) = RHS_mms%paras(:)

      ! select the T-contractor to be stored/called via C code;
      CALL fmm_select_T_con(scheme)
      ! set T_contractor pointers   
      CALL fmm_set_T_con_ptrs(Vff,RHS_mms%qlm_T)
      CALL fmm_init_T_contractors(scheme)
      CALL fmm_init_T_pair_builder(scheme,pair_type)

      ! iterate with 3D periodic repetition over supercells
      ! 4 levels equivalent = 1 + 3 + 3**2 + 3**3 + 3**4 = 121
      ! 5 levels equivalent = 364
      ncell%lo = -1093
      ncell%hi = 1093
      DO i = ncell%lo, ncell%hi 
print *, i, ncell%hi
         DO j = ncell%lo, ncell%hi 
            DO k = ncell%lo, ncell%hi 
               ! skip NN supercells
               IF ( (ABS(i)<2) .AND. (ABS(j)<2) .AND. (ABS(k)<2) )  CYCLE
               ! generate image multipoles
               DO p = 1, SIZE(RHS%raw_paras)
                  CALL pfmm_get_periodic_image(i,j,k, cell_size, cntr)
                  RHS%raw_paras(p)%cntr(:) = ref_paras(p)%cntr(:)+cntr(:)
               END DO
               ! contract with interaction tensor
               CALL fmm_gen_nonlocal_T_pairs(LHS,RHS,pair_type)
            END DO
         END DO
      END DO
      CALL fmm_close_T_pair_builder(scheme)
      CALL fmm_free_T_contractors

      CALL CPU_TIME(T1)
      CALL fmm_TIMTXT('>>> TIME USED in pfmm_get_FQ_Vff', T1-T0, IO1)
      DEALLOCATE(ref_paras)

   END SUBROUTINE pfmm_get_FQ_Vff

!-------------------------------------------------------------------------------

   SUBROUTINE pfmm_get_BQ_Vff(scheme,cell_size,LHS_raw_paras,Vff)

      IMPLICIT NONE
      TYPE(fmm_scheme),           INTENT(IN)    :: scheme
      REAL(REALK),                INTENT(IN)    :: cell_size(3)
      TYPE(raw_mm_paras), TARGET, INTENT(IN)    :: LHS_raw_paras(:)
      REAL(REALK),        TARGET, INTENT(INOUT) :: Vff(:,:)

      TYPE(gen_mm_paras)       :: LHS, RHS
      REAL(REALK), POINTER     :: qlm_T(:,:)
      REAL(REALK), POINTER     :: boxed_Vff(:,:)
      INTEGER(INTK), PARAMETER :: pair_type = LHS_box_RHS_box
      INTEGER(INTK)            :: LMAX, lm_dim
      TYPE(fmm_range) :: ncell
      REAL(REALK) :: T0, T1, cntr(3)
      INTEGER(INTK) :: i,j,k,m

      CALL CPU_TIME(T0)
      NULLIFY(qlm_T,boxed_Vff)
      NULLIFY(LHS%raw_paras,LHS%box_paras)
      NULLIFY(RHS%raw_paras,RHS%box_paras)

      CALL fmm_init_W_contractors(scheme%trans_LMAX)
      LMAX = scheme%trans_LMAX
      lm_dim = (1+LMAX)**2

      ! set up LHS and RHS boxed parameters
      CALL pfmm_get_box_paras_at_level(1,scheme,LHS%box_paras,'LHS')
      CALL pfmm_get_box_paras_at_level(1,scheme,RHS%box_paras,'RHS')
      CALL pfmm_get_box_qlm_at_level(1,scheme,qlm_T,'RHS','free')

      IF (SIZE(LHS%box_paras) > 1) CALL fmm_quit('pfmm LHS boxes!')
      IF (SIZE(RHS%box_paras) > 1) CALL fmm_quit('pfmm RHS boxes!')
      IF ( DOT_PRODUCT(RHS%box_paras(1)%cntr, RHS%box_paras(1)%cntr) > 0 ) &
         CALL fmm_quit('pfmm centre of original supercell!')

      ! allocate temporary boxed potential
      ALLOCATE (boxed_Vff(lm_dim,1))
      boxed_Vff = zero

      ! select the T-contractor to be stored/called via C code
      CALL fmm_select_T_con(scheme)
      ! set T_contractor pointers   
      CALL fmm_set_T_con_ptrs(boxed_Vff,qlm_T)
      CALL fmm_init_T_contractors(scheme)
      CALL fmm_init_T_pair_builder(scheme,pair_type)

      ! iterate with 3D periodic repetition over supercells
      m=0
      ncell%lo = -13
      ncell%hi = 13
      DO i = ncell%lo, ncell%hi 
print *, i, ncell%hi
         DO j = ncell%lo, ncell%hi 
            DO k = ncell%lo, ncell%hi 
               ! skip NN supercells
               IF ( (ABS(i)<2) .AND. (ABS(j)<2) .AND. (ABS(k)<2) ) CYCLE
               ! generate image multipoles
               CALL pfmm_get_periodic_image(i,j,k, cell_size, cntr)
               RHS%box_paras(1)%cntr(:) = cntr(:)
               ! contract with interaction tensor
               m=m+1
               CALL fmm_gen_nonlocal_T_pairs(LHS,RHS,pair_type)
            END DO
         END DO
      END DO
      print *, "# interactions = ", m
      CALL fmm_close_T_pair_builder(scheme)
      CALL fmm_free_T_contractors

      ! translate boxed potential to raw LHS centres
      CALL fmm_get_raw_Vff_from_boxed_Vff(LHS_raw_paras,scheme,boxed_Vff,Vff)
      DEALLOCATE(boxed_Vff)

      CALL CPU_TIME(T1)
      CALL fmm_TIMTXT('>>> TIME USED in pfmm_get_BQ_Vff', T1-T0, IO1)

      CALL fmm_free_W_contractors

   END SUBROUTINE pfmm_get_BQ_Vff

!-------------------------------------------------------------------------------

   SUBROUTINE pfmm_get_FMM_Vff(scheme,cell_size,LHS_raw_paras,Vff)

      IMPLICIT NONE
      TYPE(fmm_scheme),           INTENT(IN)    :: scheme
      REAL(REALK),                INTENT(IN)    :: cell_size(3)
      TYPE(raw_mm_paras), TARGET, INTENT(IN)    :: LHS_raw_paras(:)
      REAL(REALK),        TARGET, INTENT(INOUT) :: Vff(:,:)

      TYPE(gen_mm_paras)       :: LHS, RHS
      REAL(REALK), POINTER     :: qlm_T(:,:)
      REAL(REALK), POINTER     :: boxed_Vff(:,:)
      INTEGER(INTK), PARAMETER :: pair_type = LHS_box_RHS_box
      INTEGER(INTK)            :: LMAX, lm_dim, level
      TYPE(fmm_range) :: ncell
      REAL(REALK) :: T0, T1, shift(3), cntr(3)
      INTEGER(INTK) :: i,j,k,m

      CALL CPU_TIME(T0)
      NULLIFY(qlm_T,boxed_Vff)
      NULLIFY(LHS%raw_paras,LHS%box_paras)
      NULLIFY(RHS%raw_paras,RHS%box_paras)

      CALL fmm_init_W_contractors(scheme%trans_LMAX)
      ! allocate boxed potential based on number of LHS boxes
      LMAX = scheme%trans_LMAX
      lm_dim = (1+LMAX)**2
      ALLOCATE (boxed_Vff(lm_dim,1))
      boxed_Vff = zero

      ! set up LHS and RHS boxed parameters and moments
      CALL pfmm_get_box_paras_at_level(1,scheme,LHS%box_paras,'LHS')
      IF (SIZE(LHS%box_paras) > 1) CALL fmm_quit('pfmm LHS boxes!')

      ! select the T-contractor to be stored/called via C code
      CALL fmm_select_T_con(scheme)
      CALL fmm_init_T_contractors(scheme)

      m=0
      DO level = 1, pfmm_max_level
         print *, "level =", level
         ! get hierarchical boxed multipole moments
         CALL pfmm_get_box_paras_at_level(level,scheme,RHS%box_paras,'RHS')
         CALL pfmm_get_box_qlm_at_level(level,scheme,qlm_T,'RHS','free')
         IF (SIZE(RHS%box_paras) > 1) CALL fmm_quit('pfmm RHS boxes!')
!print *, "SIZE(qlm_T,1), SIZE(qlm_T,2)", SIZE(qlm_T,1), SIZE(qlm_T,2)
k=lm_dim
CALL util_matrix_norm("qlm_T norm     ",qlm_T(:k,1),k)
print *, "Total charge of aggregate supercell = ", qlm_T(1,1)
print *, "L=1 moments = "
print '(3E15.4)', qlm_T(2,1), qlm_T(3,1), qlm_T(4,1)
print *, "L=2 moments = "
print '(5E15.4)', qlm_T(5,1), qlm_T(6,1), qlm_T(7,1), qlm_T(8,1), qlm_T(9,1)
   
         ! set T_contractor pointers   
         CALL fmm_set_T_con_ptrs(boxed_Vff,qlm_T)
         CALL fmm_init_T_pair_builder(scheme,pair_type)
   
         ! iterate with 3D periodic repetition over supercells
         shift(:) = cell_size(:) * (3**(level-1))
         ncell%lo = -4
         ncell%hi = 4
         DO i = ncell%lo, ncell%hi 
            DO j = ncell%lo, ncell%hi 
               DO k = ncell%lo, ncell%hi 
                  ! skip NN supercells
                  IF ( (ABS(i)<2) .AND. (ABS(j)<2) .AND. (ABS(k)<2) ) CYCLE
                  ! generate image multipoles
                  CALL pfmm_get_periodic_image(i,j,k, shift, cntr)
                  RHS%box_paras(1)%cntr(:) = cntr(:)
                  ! contract with interaction tensor
                  m = m+1
                  CALL fmm_gen_nonlocal_T_pairs(LHS,RHS,pair_type)
               END DO
            END DO
         END DO
         CALL fmm_close_T_pair_builder(scheme)
      END DO
      CALL fmm_free_T_contractors
      print *, "# interactions = ", m

      ! translate boxed potential to raw LHS centres
      CALL fmm_get_raw_Vff_from_boxed_Vff(LHS_raw_paras,scheme,boxed_Vff,Vff)
      DEALLOCATE(boxed_Vff)

      CALL CPU_TIME(T1)
      CALL fmm_TIMTXT('>>> TIME USED in pfmm_get_FMM_Vff', T1-T0, IO1)

      CALL fmm_free_W_contractors

   END SUBROUTINE pfmm_get_FMM_Vff

!-------------------------------------------------------------------------------

END MODULE pfmm_Vff_driver

